TAUCH Lukas IF 1 TD 5
Partie 1 Data setting and Unit root test
1 data settings
Load the data
library(tseries)
library(forecast)
library(TSstudio)
library(urca)
library(tidyquant)
library(magrittr)
library(MASS)
library(eFRED)Rsxf <- fred("RSXFSN", all = FALSE)2 Modification du data et détermination du type
Rsxf_ts = ts(Rsxf$RSXFSN, start = c(1992,01), frequency = 12)
ts_plot(Rsxf_ts, title = "Sales over time" )On observe une trend avec une saisonalité qui se réperte et qui est additive car l’ampleur de la saisonnalité des fluctuations ne varient pas avec le niveau des séries chronologiques .
plot(decompose(Rsxf_ts))On observe bien une trend et la sésonalité. On ne peut donc pas appliquer un ARMA(p,q) car on observe encore une non stationnarité et une sésonalité.
acf(Rsxf_ts, lag = 120)pacf(Rsxf_ts, lag = 100 ) Notre Acf nous montre une grosse auto-corrélation donc dépendance des valeurs dans le passé ce qui montre une présence d’un AR On observe bien un ordre d’AR On observe aussi une saisonalité dans l’acf.
3 Filtrage saisonnalité
filter.outliers = function(timeseries)
{
library(foreach)
library(pracma)
## on décompose notre time series
stl_decomp = stl(timeseries, s.window = 'periodic', t.window = 13, robust = T)
## On remove la saisonnalité
TswithoutSais = seasadj(stl_decomp)
## apply a hampel filter
TsWS = hampel(TswithoutSais, k = 12)$y
## put back the seasonal component
return(TsWS)
}
Rsxf_Wseasonal <- filter.outliers(Rsxf_ts)
plot(Rsxf_Wseasonal)On a enlevé la plus part de la saisonnalité.
#Autre méthode pour enlever la saisonnalité
Rsxf_WithoutSesonal2 <- Rsxf_ts - decompose(Rsxf_ts)$seasonal
ts_plot(Rsxf_WithoutSesonal2)plot(decompose(Rsxf_WithoutSesonal2))Le même résultat
4 Check de notre Série
plot(decompose(Rsxf_Wseasonal)) On ne peut toujours pas utiliser notre signal car il reste encore la stationnarité même si on a enlevé la saisonnalité. Donc on ne peut toujours pas utilisé un modèle ARMA(p,q)
2 Unit Root tests
1 Etude stationnarité ADF
adf.test(Rsxf_Wseasonal)## Warning in adf.test(Rsxf_Wseasonal): p-value greater than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: Rsxf_Wseasonal
## Dickey-Fuller = 0.73695, Lag order = 7, p-value = 0.99
## alternative hypothesis: stationary
pvalue > O.O5 donc on ne rejette pas l’hypotese Ho qu’il y a présence de racine unitaire donc pi = 0 car (Pi = phi - 1) donc les coefficients sont significatif et Pi = 0 donc DS.
x <- ur.df(Rsxf_Wseasonal, type = "trend", lag = 2)
summary(x)##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -58234 -5318 -37 5135 76858
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3501.96094 3483.86121 1.005 0.315
## z.lag.1 -0.02160 0.02375 -0.910 0.364
## tt 29.11833 23.56172 1.236 0.217
## z.diff.lag1 -0.35715 0.05259 -6.791 4.58e-11 ***
## z.diff.lag2 -0.29702 0.05140 -5.779 1.62e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11220 on 362 degrees of freedom
## Multiple R-squared: 0.1729, Adjusted R-squared: 0.1638
## F-statistic: 18.92 on 4 and 362 DF, p-value: 3.835e-14
##
##
## Value of test-statistic is: -0.9096 5.0048 1.5297
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
#https://stats.stackexchange.com/questions/24072/interpreting-rs-ur-df-dickey-fuller-unit-root-test-resultsOn a pris lag 2 car il restait des * pour le lag = 1. Ici on voit que -0.9096 > -3.42 donc on ne rejette pas HO donc présence de racine unitaire donc \(\pi = 0\) donc la série n’est pas stationnaire. On ne rejette pas Ho quand Vtest > tau3 Ensuite 5.0048 > 4.71 ici on test Ho \(\beta1 = \beta 2 = \pi = 0\) on rejette donc Ho. On ne rejette pas Ho quand Vtest < phi2 Et enfin 1.5297 < 6.30 avec Ho \(\beta 2 = \pi = 0\) donc on ne rejette pas Ho. On ne rejette pas Ho quand Vtest < phi3
Donc la série n’est pas stationnaire avec $ = = 0 $ avec 5% d’erreur.
2 Degres d’intégration
D’après nos résultat, notre degrès d’intégration est de 2. Car le lag 1 et lag 2 sont significatifs pour le test ADF.
3 Phillips pheron
summary(ur.pp(Rsxf_Wseasonal))##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept
##
##
## Call:
## lm(formula = y ~ y.l1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48064 -7039 -404 6371 92460
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.209e+03 2.018e+03 0.599 0.549
## y.l1 1.000e+00 6.019e-03 166.131 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12320 on 367 degrees of freedom
## Multiple R-squared: 0.9869, Adjusted R-squared: 0.9868
## F-statistic: 2.76e+04 on 1 and 367 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic, type: Z-alpha is: 1.6289
##
## aux. Z statistics
## Z-tau-mu -0.1945
Ici pareil, on a une 1.6289 > -0.1945 donc on ne rejette pas H0 où \(\pi = 0\) et donc \(\phi =1\) présence de la racine unitaire et donc la série n’est pas stationnaire. On retrouve la même interprétation que Dickey et Fuller.
summary(ur.pp(diff(Rsxf_Wseasonal)))##
## ##################################
## # Phillips-Perron Unit Root Test #
## ##################################
##
## Test regression with intercept
##
##
## Call:
## lm(formula = y ~ y.l1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53527 -6369 -208 5529 85613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1523.64961 619.45790 2.460 0.0144 *
## y.l1 -0.28331 0.05019 -5.644 3.33e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11830 on 366 degrees of freedom
## Multiple R-squared: 0.08007, Adjusted R-squared: 0.07756
## F-statistic: 31.86 on 1 and 366 DF, p-value: 3.328e-08
##
##
## Value of test-statistic, type: Z-alpha is: -381.4512
##
## aux. Z statistics
## Z-tau-mu 2.9123
Ici dans notre série filtrée de la trend et des saisonnalités, on observe donc bien que -381.4512 < 2.9123 donc on rejette H0 donc la série est stationnaire.
4 KPSS
z <- ur.kpss(Rsxf_Wseasonal, type = "tau")
summary(z)##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 5 lags.
##
## Value of test-statistic is: 0.5403
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
Ici le test fait l’inverse donc 0.5403 > 0.146 donc on rejette l’hypothèse Ho donc la série n’est pas stationnaire. La série est donc bien non-stastionnaire.
5 Degres d’intégration avec KPSS
for (i in 1:7)
{
x <- summary(ur.kpss(Rsxf_Wseasonal, type = "tau", use.lag = i))
print(x)
}##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 1 lags.
##
## Value of test-statistic is: 1.4605
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 2 lags.
##
## Value of test-statistic is: 1.0109
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 3 lags.
##
## Value of test-statistic is: 0.7779
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 4 lags.
##
## Value of test-statistic is: 0.6362
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 5 lags.
##
## Value of test-statistic is: 0.5403
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 6 lags.
##
## Value of test-statistic is: 0.4716
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
##
##
## #######################
## # KPSS Unit Root Test #
## #######################
##
## Test is of type: tau with 7 lags.
##
## Value of test-statistic is: 0.4198
##
## Critical value for a significance level of:
## 10pct 5pct 2.5pct 1pct
## critical values 0.119 0.146 0.176 0.216
Avec tous les test, on observe que \(\pi = 0\) et que \(\beta1\) est significatif. Cela correspond a une non-stationnarité processus DS Et que le degré d’intégration peut-être théoriquement de 2, or on voit que différentier avec un degré d’intégration de 1 rend notre série stationnaire.
6 on applique donc le filtre avec différention 1
Rsxf_Perfect <- diff(Rsxf_Wseasonal,lag = 1, differences = 1)
ts_plot(Rsxf_Perfect)On ne voit plus de trend.
summary(ur.df(Rsxf_Perfect, type = "trend", lag = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57706 -5229 181 5063 76427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 529.76076 1178.84973 0.449 0.653
## z.lag.1 -1.68119 0.07973 -21.086 < 2e-16 ***
## tt 8.28801 5.54645 1.494 0.136
## z.diff.lag 0.30842 0.04983 6.189 1.63e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11220 on 363 degrees of freedom
## Multiple R-squared: 0.6785, Adjusted R-squared: 0.6758
## F-statistic: 255.3 on 3 and 363 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -21.0859 148.2183 222.3201
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
Ici -21.0859 << -3.42 donc rejet de Ho donc série stationnaire et donc \(\pi < 0\) Puis 148.2183 >> 4.71 donc non rejet de HO soit \(\beta 1 = \pi = 0\) Et enfin 222.3201 >> 6.30 donc non rejet de H0 soit \(\beta 1 = \beta 2 = \pi = 0\)
donc nous avons ici une série stationnaire avec un d = 1.
3 Modeliser la série
1 Most Revelant ARMA(p,q)
acf(Rsxf_Perfect, lag = 120)On voit une petite autocorrélation avec les évènements passés donc peut-être présence d’un AR.
pacf(Rsxf_Perfect, lag = 120)Avec la significativité différente de 0, on voit une corrélation différente de 0 pour potentiellement un AR(2)
On va tester le meilleur modèle par AIC et BIC :
best.order=c(0,0,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5)
{
fit.AIC=AIC(arima(Rsxf_Perfect,order=c(i,0,j)))
if(fit.AIC < best.AIC)
{
best.order=c(i,0,j)
best.arma=arima(Rsxf_Perfect,order=best.order)
best.AIC=fit.AIC
}
}
best.order ## [1] 4 0 4
best.AIC## [1] 7894.76
best.order=c(0,0,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{
fit.BIC=BIC(arima(Rsxf_Perfect,order=c(i,0,j)))
if(fit.BIC < best.BIC)
{
best.order=c(i,0,j)
best.arma=arima(Rsxf_Perfect,order=best.order)
best.BIC=fit.BIC
}
}
best.order## [1] 4 0 4
best.BIC## [1] 7933.868
On obtient pour les 2 méthodes un ARMA(4,4).
Pour une alternative, on peut modéliser notre méthode avant le filtrage pour la rendre stationnaire avec un ARIMA(p,d,q), on va donc faire d’une pierre deux coups. Rendre la série stationnaire et lui trouver un bon model.
2 Test du model
Estimation <- arima(Rsxf_Perfect,order=c(4,0,4))
values_observed <- Rsxf_Perfect
residus <- Estimation$residuals
values_estimated <-(values_observed - residus)par(mfrow=c(2,1))
plot(values_estimated,type='l',ylab="values_estimated")
plot(Rsxf_Perfect, type = "l")
On voit une ressemblance même s’il le modèle estimé ne fit pas
parfaitement.
Quality Checks : Residus
Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:20)
{
Quality_Checks[i,1] = Box.test(residus,lag=i,type="Ljung")$statistic
Quality_Checks[i,2] = Box.test(residus,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])plot(Quality_Checks[,2])On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle. Cependant, Nous avons donc manqué quelque chose dans la modélisation de notre model. Peut-être le fait d’avoir fait une différence avec un degré d’uintégration à 1 au lieux de 2.
acf(residus)pacf(residus)On retrouve dans notre ACF, une petite autocorrélation mais ce modèle fit le mieux avec notre data.
qqnorm(residus)
qqline(residus)On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.
plot(density(residus))On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.
QUALITY CHECKS Significativité des coefficients trouvés
\[RMSE = \sqrt{\frac{\sum_{i = 1}^{N} (Predicted_{i} - Actual_{i})^{2}}{N} }\]
RMSE <- function(x,y)
{
# with x is the predicted i vector
# y the actual i vector
rmse <- sqrt(mean((x-y)^2))
return(rmse)
}RMSE(values_estimated,values_observed)## [1] 10395.05
Exercie 4 Estimating ARIMA(p,d,q)
jnj = tq_get("JNJ", get="stock.prices", from ="1997-01-01") %>%
tq_transmute(mutate_fun = to.period, period = "months")head(jnj)## # A tibble: 6 × 7
## date open high low close volume adjusted
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1997-01-31 28.6 29.2 28.5 28.9 5477800 15.5
## 2 1997-02-28 28.6 29.2 28.3 28.8 4997200 15.4
## 3 1997-03-31 27.5 27.7 26.4 26.4 8000400 14.2
## 4 1997-04-30 30 30.8 30 30.6 4504600 16.4
## 5 1997-05-30 29.4 30.1 29.2 30 4558800 16.2
## 6 1997-06-30 32.4 32.4 31.6 32.2 5272000 17.3
ts_plot(jnj)acf(jnj$close, lag = 100)pacf(jnj$close, lag = 100)on observe une trend très important au niveau du close. J’étudie le close pour le Johnson & Johnson stock price car la valeur du close représente bien une date précise. De plus le close, représente tout aussi bien le low et high un peu comme une moyenne sur une journée donnée. Ici, on voit que sur l’acf qu’il y a présence d’un AR car AR = MA(Infini)
1 Détermination degré d’intégration
summary(ur.df(jnj$close, type = "trend", lag = 2))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.3341 -2.0488 -0.0715 2.3043 14.5427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.561670 0.560610 1.002 0.317193
## z.lag.1 -0.021627 0.015865 -1.363 0.173828
## tt 0.011900 0.007037 1.691 0.091858 .
## z.diff.lag1 -0.094066 0.056652 -1.660 0.097861 .
## z.diff.lag2 -0.208921 0.056455 -3.701 0.000255 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.163 on 304 degrees of freedom
## Multiple R-squared: 0.06368, Adjusted R-squared: 0.05136
## F-statistic: 5.169 on 4 and 304 DF, p-value: 0.0004839
##
##
## Value of test-statistic is: -1.3632 3.3917 1.5724
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
On observe bien que -1,3 > -3,42 donc on ne rejette pas H0 soit \(\pi = 0\) donc série non stationnaire. J’ai fais la l’ordre 2 car on a une p-value tres significative 3* avec une degré 2. On suppose donc un degré d’intégration de 2.
summary(ur.df(diff(jnj$close), type = "trend", lag = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.6654 -2.0530 -0.0332 2.2828 14.3702
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.168653 0.477828 0.353 0.7244
## z.lag.1 -1.327957 0.082388 -16.118 <2e-16 ***
## tt 0.003020 0.002666 1.133 0.2582
## z.diff.lag 0.220373 0.055905 3.942 0.0001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.169 on 305 degrees of freedom
## Multiple R-squared: 0.5663, Adjusted R-squared: 0.562
## F-statistic: 132.7 on 3 and 305 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -16.1183 86.6028 129.9037
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
Or ici j’intègre avec un degré de 1 et cela me renvoie a une série stationnaire car -16 < -3 donc on rejette H0 => donc série stationnaire.
Au final le degré d’intégration est de 1.
2 Order ARIMA(p,d,q)
On connait d = 1 donc on peut faire la méthode du AIC et BIC pour déterminer l’ordre du ARIMA.
best.order=c(0,1,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5)
{
fit.AIC=AIC(arima(jnj$close,order=c(i,1,j)))
if(fit.AIC < best.AIC)
{
best.order=c(i,1,j)
best.arma=arima(jnj$close,order=best.order)
best.AIC=fit.AIC
}
}
best.order ## [1] 4 1 5
best.AIC## [1] 1772.018
best.order=c(0,1,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{
fit.BIC=BIC(arima(jnj$close,order=c(i,1,j)))
if(fit.BIC < best.BIC)
{
best.order=c(i,1,j)
best.arma=arima(jnj$close,order=best.order)
best.BIC=fit.BIC
}
}
best.order## [1] 0 1 2
best.BIC## [1] 1789.379
On se retrouve avec 2 modeles possible soit : AIC ARIMA(5,1,5) BIC ARIMA(0,1,2)
Cependant, avec notre PACF nous avions vu qu’il n’y avait pas d’ordre pour les AR donc le BIC semble être le plus représentatif.
auto.arima(jnj$close)## Series: jnj$close
## ARIMA(2,1,2) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 drift
## -0.2286 0.4235 0.0977 -0.6762 0.4724
## s.e. 0.1591 0.1557 0.1329 0.1322 0.1225
##
## sigma^2 = 17.05: log likelihood = -879.94
## AIC=1771.87 AICc=1772.15 BIC=1794.31
Avec une auutre méthode, on retrouve le modèle ARIMA(0,1,2) comme étant le plus significatif.
3 Fit ARIMA estimated valued
Estimation <- arima(jnj$close,order=c(0,1,2))
values_observed <- jnj$close
residus <- Estimation$residuals
values_estimated <-(values_observed - residus)par(mfrow=c(2,1))
plot(values_estimated,type='l',ylab="values_estimated")
plot(jnj$close, type = "l")On observe une grosse ressemblence.
4 QUality checks on Residuals
Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
Quality_Checks[i,1] = Box.test(residus,lag=i,type="Ljung")$statistic
Quality_Checks[i,2] = Box.test(residus,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])plot(Quality_Checks[,2])On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle. Cependant, on a un drop de pvalue donc nous avons donc manqué quelque chose dans la modélisation de notre model. Peut-être le fait d’avoir fait une différence avec un degré d’uintégration à 1 au lieux de 2.
acf(residus)#pacf(residus)On voit dans l’ACF de nos résidus sont indépendants les uns des autres.
qqnorm(residus)
qqline(residus)On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.
plot(density(residus))On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.
5 Forecast
predic <-forecast(values_estimated,level=c(95),h=10)
plot(predic, xlim=c(250,320))
prediction<- predict(predic,jnj$close,se = TRUE,interval="prediction")
lines(prediction$lower,lty="dashed",col="red")
lines(prediction$upper,lty="dashed",col="red")CALULER CONFIDENCE DU FORCAST
5 Unit root test
1 Zivot and Andrews
Le test de Zivot et Andrews améliore le test de stationnarité d’une série temporelle. En effet, cette méthode dait un test de racine unitaire avec une rupture structurelle. Cette rupture est estimé.
Soit notre Hypothèse H0 : \(\pi = 0\) donc série non stationnaire H1 : \(\pi < 0\) avec soit : - 0 break - 1 break -> constante -> Trend -> constante + Trend
2 Generer 3 marches aléatoires
set.seed(123)
ut <- rnorm(500,0,1)
t <- 1:500
# Random walk
y <- cumsum(ut)
# Random walk with break
ybreak <- c()
for (i in 1:500)
{
ifelse( i< 250, ybreak[i] <- cumsum(ut[i]),ybreak[i] <- 20 + cumsum(ut[i]))
}
#random walk with break and trend
ytrend <- c()
for (i in 1:500)
{
ifelse( i< 250, ytrend[i] <- cumsum(ut[i]) + 0.05*t[i] ,ytrend[i] <- 20 + cumsum(ut[i]) + 0.05*t[i])
}
plot(ytrend, type = "l")
lines(ybreak, type = "l", col = "blue")
lines(y, type = "l", col = "red")#library(strucchange)
#library(zoo)
#nbreak <- breakpoints(ybreak~t, h = 20) 3 Test Zivot and Andrews
Dans un premier temps on va regarder avec un test adf avec comme condition On ne rejette pas Ho quand Vtest > tau3 donc série non stationnaire.
#summary(ur.df(y, type = "trend", lag = 1))
#summary(ur.df(ybreak, type = "trend", lag = 1))
#summary(ur.df(ytrend, type = "trend", lag = 1))
# J'évite de le print mais je vous met une interpretation en dessous.Pour y la marche aléatoire normale : -2.41 > -3.42 donc la série n’est pas stationnaire De même pôur la marche aléatoire avec un break (constante) : -2.6139 > -3.42 Et pour la marche aléatoire avec break et trend -2.6139 > -3.42
Donc d’après notre test ADF aucune des 3 n’est stationnaires.
Maintennant on applique le test de Zivot et Andrews.
summary(ur.za(y, model = "both", lag = 1))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.95301 -0.59372 -0.02188 0.61290 2.99265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2298772 0.1288764 1.784 0.07509 .
## y.l1 0.9486190 0.0140338 67.595 < 2e-16 ***
## trend -0.0012935 0.0008162 -1.585 0.11366
## y.dl1 -0.0347155 0.0448744 -0.774 0.43953
## du 0.4979822 0.1917264 2.597 0.00968 **
## dt 0.0023018 0.0013215 1.742 0.08217 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.963 on 492 degrees of freedom
## (2 observations effacées parce que manquantes)
## Multiple R-squared: 0.9613, Adjusted R-squared: 0.9609
## F-statistic: 2445 on 5 and 492 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -3.6612
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 263
summary(ur.za(ybreak, model = "both", lag = 1))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7389 -0.6166 -0.0155 0.6387 3.2192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0736376 0.1261368 0.584 0.560
## y.l1 -0.0444094 0.0363614 -1.221 0.223
## trend -0.0006451 0.0008713 -0.740 0.459
## y.dl1 0.0092039 0.0287074 0.321 0.749
## du 21.0483874 0.7376339 28.535 <2e-16 ***
## dt 0.0006860 0.0012232 0.561 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9762 on 492 degrees of freedom
## (2 observations effacées parce que manquantes)
## Multiple R-squared: 0.9908, Adjusted R-squared: 0.9907
## F-statistic: 1.054e+04 on 5 and 492 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -28.723
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 249
summary(ur.za(ytrend, model = "both", lag = 1))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7389 -0.6166 -0.0155 0.6387 3.2192
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.070957 0.126134 0.563 0.574
## y.l1 -0.044409 0.036361 -1.221 0.223
## trend 0.051575 0.002001 25.780 <2e-16 ***
## y.dl1 0.009204 0.028707 0.321 0.749
## du 21.048387 0.737634 28.535 <2e-16 ***
## dt 0.000686 0.001223 0.561 0.575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9762 on 492 degrees of freedom
## (2 observations effacées parce que manquantes)
## Multiple R-squared: 0.9966, Adjusted R-squared: 0.9966
## F-statistic: 2.899e+04 on 5 and 492 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -28.723
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 249
Pour la marche aléatoire simple : -3.6612 > -5.08 donc non rejet de H0 donc \(\pi = 0\) donc non stationaire avec potentiellement un break à l’indice 263.
Cependant pour la marche aléatoire avec break (en constante) : -28.723 < -5.08 donc on rejette H0 donc la série est stationnaire avec un break à l’indice 249, pour rappelle j’ai mis mon break à l’indice 250. Le test est donc assez précis.
Et enfin pour la marche aléatoire avec break et trend, on a -28.723 < -5.08 donc de même série stationnaire mais avec un break à l’indice 249.
on voit donc l’importance de faire des tests avec ruptures structurelles.
4 ZA in US retail sales
on replot notre jeu de donnée sans séonalité.
ts_plot(Rsxf_Wseasonal)On observe potentiellement un break au alentour des années 2007.
#summary(ur.df(Rsxf_Wseasonal, type = "trend", lag = 1))
summary(ur.za(Rsxf_Wseasonal, model = "both", lag = 1))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57499 -6186 314 5521 54650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.251e+04 6.053e+03 7.023 1.08e-11 ***
## y.l1 7.176e-01 4.024e-02 17.833 < 2e-16 ***
## trend 2.515e+02 3.605e+01 6.977 1.44e-11 ***
## y.dl1 -1.607e-01 5.096e-02 -3.153 0.00175 **
## du 3.749e+04 5.807e+03 6.457 3.45e-10 ***
## dt -2.886e+02 4.486e+02 -0.643 0.52048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10940 on 362 degrees of freedom
## (2 observations effacées parce que manquantes)
## Multiple R-squared: 0.9897, Adjusted R-squared: 0.9896
## F-statistic: 6972 on 5 and 362 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -7.0187
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 350
Je rappelle qu’avec DF on trouve : -2.3639 > -3.42 donc séries non stationnaire Alors qu’avec ZA on trouve : -7.0187 < -5.08 donc série stationnaire avec un break à l’indice 350 Donc on peut se poser la question si ce break, chute en 2007, à l’indice 350 ne faussait pas tous nos tests.
6 Modeling the business cycle
Load data
Gdp <- fred("GDP", all = FALSE)
Gdp <- na.omit(Gdp, "GDP")J’enlève toute valeure null.
Modification en TS
Gdp_ts = ts(Gdp$GDP, start = c(1947,01), frequency = 4)
ts_plot(Gdp_ts, title = "Gross Domestic Product over time" )On observe un belle trend avec pas de saisonalité apparente.
Tests sur notre Time serie
#plot(decompose(Gdp_ts))acf(Gdp_ts, lag = 100)pacf(Gdp_ts)ACF -> Grosse autocorrélation, peut-être présence d’un AR. PACF -> 1 seul corrélation à 1
summary(ur.df(Gdp_ts, type = "trend", lag = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2190.79 -10.84 -1.27 19.13 1231.84
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.007254 26.986428 0.000 0.9998
## z.lag.1 0.009917 0.004156 2.386 0.0176 *
## tt 0.184980 0.323210 0.572 0.5675
## z.diff.lag -0.128714 0.058249 -2.210 0.0279 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 170.4 on 297 degrees of freedom
## Multiple R-squared: 0.1706, Adjusted R-squared: 0.1623
## F-statistic: 20.37 on 3 and 297 DF, p-value: 4.955e-12
##
##
## Value of test-statistic is: 2.3865 38.4439 30.0376
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
On observe avec le test de Dickey Fuller que 2.3865 > -3.42 donc non rejet de H0 -> \(\pi = 0\) donc série non stationnaire. Par rapport au drift et trend notre T-stat sont supérieures aux valeurs 5% donc on rejette leur H0 donc rejet de nullité.
Mais nous aller voir si ce n’est pas un problème de break
summary(ur.za(Gdp_ts, model = "both", lag = 1))##
## ################################
## # Zivot-Andrews Unit Root Test #
## ################################
##
##
## Call:
## lm(formula = testmat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2078.14 -21.93 0.01 30.15 380.78
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -34.196952 23.047319 -1.484 0.138937
## y.l1 0.995171 0.003762 264.545 < 2e-16 ***
## trend 0.969200 0.282956 3.425 0.000701 ***
## y.dl1 -0.178682 0.051268 -3.485 0.000566 ***
## du 909.832281 106.208351 8.566 6.03e-16 ***
## dt -61.533073 19.202554 -3.204 0.001502 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 141.9 on 295 degrees of freedom
## (2 observations effacées parce que manquantes)
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 1.439e+05 on 5 and 295 DF, p-value: < 2.2e-16
##
##
## Teststatistic: -1.2838
## Critical values: 0.01= -5.57 0.05= -5.08 0.1= -4.82
##
## Potential break point at position: 294
On a -1.2838 > -5.08 donc non rejet de H0 qui comme adf \(\pi = 0\) donc série non sationnaire mais potentielle break à l’indice 294.
Conclusion : On sait donc que notre série est bien non sattionnaire avec d = 2
Différentiation et estimation du model ARIMA
Différention avec d = 1
Gdp_perfect <- diff(Gdp_ts, differences = 2)
ts_plot(Gdp_perfect)Test de stationnarité
summary(ur.df(Gdp_perfect, type = "trend", lag = 1))##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2055.74 -18.99 1.60 14.27 2067.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.37458 22.33321 -0.151 0.880
## z.lag.1 -2.19791 0.09628 -22.828 < 2e-16 ***
## tt 0.04451 0.12842 0.347 0.729
## z.diff.lag 0.37815 0.05392 7.013 1.59e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 191.6 on 295 degrees of freedom
## Multiple R-squared: 0.8264, Adjusted R-squared: 0.8246
## F-statistic: 468 on 3 and 295 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -22.8284 173.7121 260.5682
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.98 -3.42 -3.13
## phi2 6.15 4.71 4.05
## phi3 8.34 6.30 5.36
#summary(ur.pp(Gdp_perfect))
#summary(ur.za(Gdp_perfect, model = "both", lag = 1))avec notre test on a -22.8284 < -3.42 donc rejet de H0 donc série stationnaire.
Estimation du model ARIMA avec d = 2 et ARMA
acf(Gdp_perfect)pacf(Gdp_perfect)ACF -> Très petit autocorrélation mais existante MA(1) PACF -> petite corrélation AR(1)
auto.arima(Gdp_ts)## Series: Gdp_ts
## ARIMA(0,2,3)
##
## Coefficients:
## ma1 ma2 ma3
## -1.0690 0.2205 -0.0932
## s.e. 0.0579 0.0859 0.0577
##
## sigma^2 = 30054: log likelihood = -1978.49
## AIC=3964.97 AICc=3965.11 BIC=3979.8
auto.arima(Gdp_perfect)## Series: Gdp_perfect
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## -1.0700 0.0990 0.9300
## s.e. 0.0533 0.0575 0.3959
##
## sigma^2 = 29847: log likelihood = -1977.81
## AIC=3963.61 AICc=3963.75 BIC=3978.44
ARIMA(p,2,q)
best.order=c(0,2,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5)
{
fit.AIC=AIC(arima(Gdp_ts,order=c(i,2,j), method = "ML"))
if(fit.AIC < best.AIC)
{
best.order=c(i,2,j)
best.arma=arima(Gdp_ts,order=best.order, method = "ML")
best.AIC=fit.AIC
}
}
best.order ## [1] 1 2 1
best.AIC## [1] 3964.754
best.order=c(0,2,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{
fit.BIC=BIC(arima(Gdp_ts,order=c(i,2,j), method = "ML"))
if(fit.BIC < best.BIC)
{
best.order=c(i,2,j)
best.arma=arima(Gdp_ts,order=best.order, method = "ML")
best.BIC=fit.BIC
}
}
best.order## [1] 0 2 1
best.BIC## [1] 3975.57
Avec le AIC, on obtient un ARIMA(1,2,1) et pour le BIC un ARIMA(0,2,1)
Determination ARMA(p,q) avec notre série différentié
best.order=c(0,0,0)
best.AIC=10000
for (i in 0:5) for (j in 0:5)
{
fit.AIC=AIC(arima(Gdp_perfect,order=c(i,0,j), method = "ML"))
if(fit.AIC < best.AIC)
{
best.order=c(i,0,j)
best.arma=arima(Gdp_perfect,order=best.order, method = "ML")
best.AIC=fit.AIC
}
}
best.order ## [1] 0 0 3
best.AIC## [1] 3961.523
best.order=c(0,0,0)
best.BIC= 10000
for (i in 0:5) for (j in 0:5)
{
fit.BIC=BIC(arima(Gdp_perfect,order=c(i,0,j), method = "ML"))
if(fit.BIC < best.BIC)
{
best.order=c(i,0,j)
best.arma=arima(Gdp_perfect,order=best.order, method = "ML")
best.BIC=fit.BIC
}
}
best.order## [1] 0 0 1
best.BIC## [1] 3975.339
Avec le AIC, on obtient un ARIMA(0,0,3) et pour le BIC un ARIMA(0,0,1)
En conclusion, j’utilise le ARIMA(1,2,1) car il est correct par rapport à notre ACF et PACF Et l’ARMA(0,0,2)
Valeurs estimées avec nos 2 modèles ARIMA(1,2,1) et ARMA(0,0,2)
ARIMA(1,2,1)
Estimation_Gdp <- arima(Gdp_ts,order=c(1,2,1), method = "ML")
values_observed_Gdp <- Gdp_ts
residus_Gdp <- Estimation_Gdp$residuals
values_estimated_Gdp <-(values_observed_Gdp - residus_Gdp)par(mfrow=c(2,1))
plot(values_estimated_Gdp,type='l',ylab="values_estimated")
plot(Gdp_ts, type = "l")ARMA(0,0,2)
Estimation_Gdp2 <- arima(Gdp_perfect,order=c(0,0,2), method = "ML")
values_observed_Gdp2 <- Gdp_perfect
residus_Gdp2 <- Estimation_Gdp2$residuals
values_estimated_Gdp2 <-(values_observed_Gdp2 - residus_Gdp2)par(mfrow=c(2,1))
plot(values_estimated_Gdp2,type='l',ylab="values_estimated")
plot(Gdp_perfect, type = "l")Quality checks
ARIMA(1,2,1)
Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
Quality_Checks[i,1] = Box.test(residus_Gdp,lag=i,type="Ljung")$statistic
Quality_Checks[i,2] = Box.test(residus_Gdp,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])plot(Quality_Checks[,2])On observe que la p-value > 0,05 jusqu’à un lag < 10 les résidus sont donc indépendant sur cette intervalle. Nous avons donc manqué quelque chose dans la modélisation de notre model.
acf(residus_Gdp)#pacf(residus_Gdp)On voit dans l’ACF de nos résidus sont indépendants les uns des autres.
qqnorm(residus_Gdp)
qqline(residus_Gdp)On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.
plot(density(residus_Gdp))On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.
RMSE(values_estimated_Gdp,values_observed_Gdp)## [1] 172.4497
ARMA(0,0,2)
Quality_Checks = matrix(0,ncol=2,nrow=20)
for(i in 1:10)
{
Quality_Checks[i,1] = Box.test(residus_Gdp2,lag=i,type="Ljung")$statistic
Quality_Checks[i,2] = Box.test(residus_Gdp2,lag=i,type="Ljung")$p.value
}
plot(Quality_Checks[,1])plot(Quality_Checks[,2])On observe que la p-value > 0,05 les résidus sont donc indépendant sur cette intervalle.
acf(residus_Gdp2)#pacf(residus_Gdp2)On voit dans l’ACF de nos résidus sont indépendants les uns des autres.
qqnorm(residus_Gdp2)
qqline(residus_Gdp2)On observe bien que les résidus suivent une distribution normale car la plupart des points se situent sur la ligne.
plot(density(residus_Gdp2))On voit bien ici que les résidus trouvés forment bien une densité normale centré en 0.
RMSE(values_estimated_Gdp2,values_observed_Gdp2)## [1] 171.9032
On remarque que les résidus des 2 modèles sont très similaires en terme de leur forme et la même RMSE
Forecasting
ARIMA(1,2,1)
predic_Gdp <-forecast(values_estimated_Gdp,level=c(95),h=10)
plot(predic_Gdp, xlim=c(2000,2025))
prediction_Gdp<- predict(predic_Gdp,Gdp_ts,se.fit = TRUE,interval="prediction")
lines(prediction_Gdp$lower,lty="dashed",col="red")
lines(prediction_Gdp$upper,lty="dashed",col="red")ARMA(0,0,2)
predic_Gdp2 <-forecast(Estimation_Gdp2,level=c(95),h=5)
plot(predic_Gdp2,xlim=c(2015,2024))
prediction_Gdp2<- predict(predic_Gdp2,Gdp_perfect,se.fit = TRUE, interval="prediction")
lines(prediction_Gdp2$lower,lty="dashed",col="red")
lines(prediction_Gdp2$upper,lty="dashed",col="red")Voici les différents forecasting. On remarque que le forecasting du model ARIMA(1,2,1) semble meilleur que celui de ARMA(0,0,2).